1. Stats de base, modèles linéaires, GLM
Distribution statistique
Décrit les probabilités d’apparition des valeurs d’une distribution donnée.
Distribution normale
On parle aussi souvent d’une distribution Gaussienne.
\[ f(x) = \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2} \]
où \(\mu\) et \(\sigma\) sont des paramètres, respectivement la moyenne et l’écart-type.
Forme de la distribution normale
x <- seq(-5, 15, by = 0.01)
mu <- 5
sigma <- 2
y <- (1/(sigma * sqrt(2 * pi)) * exp((-1/2) * ((x - mu)/sigma)^2))
plot(x, y, type = "l")
Simuler une distribution normale
On peut utiliser la fonction rnorm pour simuler des
valeurs provenant d’une distribution normale.
set.seed(123)
n <- 1000
x <- rnorm(n, mu, sigma)
hist(x, breaks = 50)
Interprétation
v <- seq(-5, 15, by = 0.01)
hist(x, breaks = 50, freq = FALSE, xlim = range(v))
curve(dnorm(x, mu, sigma), add = TRUE)
y <- dnorm(v, mu, sigma)
polygon(c(v[v >= 6 & v <= 8], 8, 6), c(y[v >= 6 & v <= 8], 0, 0), col = alpha("red", 0.5), border = NA)
L’aire sous la courbe est de 1 et la proportion de la surface en
rouge représente la probabilité d’obtenir une valeur entre 6 et 8 avec
une moyenne de \(\mu = 5\) et un
écart-type de \(\sigma = 2\). On peut
estimer cette probabilité avec la fonction pnorm intégrée à
R. Ce type de fonction est intégré pour plusieurs distribution
statistique (voir ?Distributions).
Probabilité d’obtenir une valeur?
pnorm(5, mean = 5, sd = 2) # probabilité d'obtenir une valeur <= 5## [1] 0.5pnorm(-2, mean = 5, sd = 2) # probabilité d'obtenir une valeur <= -2## [1] 0.0002326291pnorm(4, mean = 5, sd = 2) # probabilité d'obtenir une valeur <= 4## [1] 0.3085375Probabilité d’obtenir une valeur entre 6 et 8
pnorm(8, 5, 2) - pnorm(6, 5, 2)## [1] 0.2417303Distribution d’une variable
\[ X \sim \mathcal{N}(\mu,\,\sigma^{2})\,\]
\(X\) suit une distribution normale avec une moyenne \(\mu\) et une variance \(\sigma^{2}\)
Erreur-type
\[ se = \frac{sd}{\sqrt{n}} \]
ou
\[ \sigma_{\bar{x}} = \frac{\sigma_{x}}{\sqrt{n}} \]
n <- 1e+05
pop <- rnorm(n, 100, 20) # population fictive de n individus avec une moyenne de 100 et écart-type de 20
ech <- sample(pop, 30) # échantillon aléatoire de 30 individus
mean(ech) # moyenne## [1] 100.5529sd(ech)/sqrt(30) # erreur-type## [1] 3.27161Simulation
Ici, on refait cet échantillonnage nreps fois et on
calcule l’écart-type des différentes moyennes obtenues.
nreps <- 1000
ech1000 <- sapply(1:nreps, function(i) {
s <- sample(pop, 30)
mean(s)
})
hist(ech1000)
Écart-type des valeurs simulées
sd(ech1000)## [1] 3.734432Comparons ceci à l’erreur-type calculé à partir de notre seul échantillon:
sd(ech)/sqrt(30) # erreur-type## [1] 3.27161Interprétation de l’erreur-type
L’écart-type de ces moyennes correspond à l’erreur-type estimée à partir d’un seul échantillon. En d’autres mots, l’erreur-type représente:
l’écart-type si des moyennes obtenues si on refaisait plusieurs fois le même échantillonnage
une mesure de précision dans l’estimation de la moyenne
Le concept d’erreur-type s’applique aussi à d’autres paramètres estimés et non seulement à la moyenne.
L’erreur-type d’une statistique ou d’un paramètre est une estimation de l’écart-type de sa distribution d’échantillonnage.
Intervalle de confiance
Reprenons notre population fictive pop qui a une moyenne
de 100 et un écart-type de 20. On se rappelle qu’un intervalle de
confiance pour une moyenne provenant d’un échantillon distribué
normalement (ou presque) se calcule de la façon suivante:
\[\bar{x} ± t_{dl,\alpha/2}\sigma_{\bar{x}}\]
On peut calculer cet intervalle de confiance avec R:
mean(ech) + qt(0.025, df = length(ech) - 1, lower.tail = FALSE) * (sd(ech)/sqrt(length(ech))) * c(-1, 1)## [1] 93.86168 107.24407Simulation
Maintenant, reprenons note population fictive et calculons à chaque fois cet intervalle de confiance. Ensuite, calculons quelle est la proportion des intervalles qui contiennent la véritable moyenne de 100.
nreps <- 1000
ech1000 <- lapply(1:nreps, function(i) {
s <- sample(pop, 30)
c(mean(s), mean(s) + qt(0.025, df = length(s) - 1, lower.tail = FALSE) * (sd(s)/sqrt(length(s))) * c(-1, 1))
})
ci <- as.data.frame(do.call("rbind", ech1000))
names(ci) <- c("mean", "lowerCI", "upperCI")
sum(100 >= ci$lowerCI & 100 <= ci$upperCI)/nreps## [1] 0.952Interprétation de l’intervalle de confiance
Un intervalle de confiance à 95% veut dire que 95% des intervalles de confiance qu’on obtiendrait si on refaisait le même échantillonnage contiendront la véritable moyenne de la population.
Plus généralement, on devrait parler d’un paramètre (moyenne, variance, pente, etc.) plutôt que d’une moyenne.
On entend souvent à tort que cela veut dire qu’il y a 95% de chances que le paramètre soit à l’intérieur des bornes. Or le paramètre est considéré comme étant fixe et la probabilité porte sur le fait que l’intervalle contienne la véritable valeur du paramètre ou pas.
Tout comme l’erreur-type, l’intervalle de confiance est une mesure de précision dans l’estimation d’un paramètre.
Valeur de p
Un exemple avec le test de \(t\) et une population fictive.
set.seed(1)
n <- 10000
pop1 <- rnorm(n, 20, 3) # population fictive de n individus avec une moyenne de 20
pop2 <- rnorm(n, 21, 3) # population fictive de n individus avec une moyenne de 21
ech1 <- sample(pop1, 10) # échantillon de 10 individus de chaque population
ech2 <- sample(pop2, 10)
test <- t.test(ech1, ech2)
test##
## Welch Two Sample t-test
##
## data: ech1 and ech2
## t = -1.189, df = 14.686, p-value = 0.2533
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -4.030060 1.147318
## sample estimates:
## mean of x mean of y
## 19.22124 20.66261Interprétation de la Valeur de p
plot(0, 0, xlim = c(-4, 4), ylim = c(0, 0.5), type = "n")
curve(dt(x, test$parameter), add = TRUE)
v <- seq(-5, 15, by = 0.01)
y <- dt(v, test$parameter)
polygon(c(v[v <= test$statistic], test$statistic, min(v)), c(y[v <= test$statistic], 0, 0), col = alpha("red", 0.5), border = NA)
polygon(c(v[v >= abs(test$statistic)], max(v), abs(test$statistic)), c(y[v >= abs(test$statistic)], 0, 0), col = alpha("red", 0.5), border = NA)
abline(v = rep(test$statistic, 2) * c(1, -1), lwd = 2, lty = 3, col = "red")
Simulations
Reprenons le test de \(t\) précédent, mais cette fois, prenons nos échantillons dans la même population (pop1) à chaque fois.
nreps <- 1000
pvalue <- sapply(1:nreps, function(i) {
ech1 <- sample(pop1, 5)
ech2 <- sample(pop1, 5)
test <- t.test(ech1, ech2)
test$p.value #statistic
})
hist(pvalue, breaks = 50)
abline(v = 0.05, col = "red", lwd = 2, lty = 3)
sum(pvalue <= 0.05)/length(pvalue)## [1] 0.051Modèles linéaires simples
On a ici la formulation classique ou habituelle d’un modèle linéaire simple.
\[ y = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{2} + ... + \beta_{n}x_{n} + \epsilon \]
\[ \epsilon \sim \mathcal{N}(0,\,\sigma^{2})\ \]
Autre formulation
\[ \mu = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{2} + ... + \beta_{n}x_{n}\]
\[ y \sim \mathcal{N}(\mu,\,\sigma^{2})\ \]
où
\(y\) = observations
\(\mu\) = prédicteur linéaire
On a donc un modèle représentant le lien entre nos variables et on cherche à estimer les paramètres de ce modèle à l’aide de données.
Suppositions de bases
Quelles sont les suppositions de bases d’un modèle linéaire tel que celui présenté?
Normalité des résidus?
Homoscédasticité de la variance (la variance ne dépend pas de ou ne varie pas avec \(x_{n}\))
Linéarité des relations entre \(x_{n}\) et \(y\)?
Corrélation entre les \(x_{n}\)?
Indépendence entre les observations?
Simulations
Une des meilleures façons de s’assurer qu’on comprend comment fonctionne un modèle
Permet d’étudier le comportement d’un modèle avec des valeurs connues
Simulons un modèle
# nb d'observations
n <- 75
# valeurs fictives des variables explicatives selon une distribution uniforme
x1 <- runif(n, 0, 100)
x2 <- runif(n, 0, 10)
x3 <- runif(n, 0, 1)
# valeurs des paramètres beta
b0 <- 100
b1 <- 0.2
b2 <- 1
b3 <- -10
# prédicteur linéaire
mu <- b0 + b1 * x1 + b2 * x2 + b3 * x3
# observations à partir d'une distribution normale avec un écart-type de 5
y <- rnorm(n, mu, 5)
# données
d <- data.frame(y, x1, x2, x3)Fonction lm
La fonction lm permet de faire des modèles linéaires
simples, ce qui inclue les régressions simples, les ANOVAs, ANCOVAs,
régression multiples, etc.
Remarquez l’utilisation de l’argument data où la
fonction ira chercher les différents éléments spécifiés dans la
formule.
m <- lm(y ~ x1 + x2 + x3, data = d)Coefficients
On peut appliquer plusieurs fonctions pour extraire les éléments d’un
modèle, en particulier un modèle effectué avec la fonction
lm.
coef(m)## (Intercept) x1 x2 x3
## 96.9666068 0.1882304 1.0717263 -4.0315717Sommaire
La fonction summary est la plus utile pour extraire
l’information d’un modèle.
summary(m)##
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.518 -2.927 0.206 2.739 12.382
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 96.96661 1.80048 53.856 < 2e-16 ***
## x1 0.18823 0.01909 9.858 6.23e-15 ***
## x2 1.07173 0.18827 5.693 2.62e-07 ***
## x3 -4.03157 1.92616 -2.093 0.0399 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.02 on 71 degrees of freedom
## Multiple R-squared: 0.6444, Adjusted R-squared: 0.6294
## F-statistic: 42.89 on 3 and 71 DF, p-value: 6.301e-16Vérification des suppositions de bases
On peut utiliser la fonction plot qui lorsqu’elle est
appliquée sur un objet de classe lm produira des graphiques
permettant de vérifier certaines suppositions de bases ou conditions
d’application.
par(mfrow = c(2, 2))
plot(m)
Fonctions accessoires
On peut également créer facilement certains de ces graphiques avec d’autres fonctions.
par(mfrow = c(1, 3))
plot(resid(m), fitted(m))
hist(resid(m))
qqnorm(resid(m))
qqline(resid(m))
Facteurs
# simulation d'un facteur
d$fac <- factor(sample(c("A", "B", "C"), nrow(d), replace = TRUE))
d$y <- d$y + 2 * as.integer(d$fac) - 1
# modèle modifié avec le facteur
m <- lm(y ~ x1 + x2 + x3 + fac, data = d)
summary(m)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + fac, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.3549 -3.0036 0.3166 2.6305 13.2288
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 97.2099 1.9870 48.924 < 2e-16 ***
## x1 0.1861 0.0190 9.794 1.09e-14 ***
## x2 1.0464 0.1854 5.645 3.39e-07 ***
## x3 -3.8057 1.8956 -2.008 0.0486 *
## facB 1.9366 1.4416 1.343 0.1836
## facC 6.5096 1.3744 4.736 1.12e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.932 on 69 degrees of freedom
## Multiple R-squared: 0.696, Adjusted R-squared: 0.674
## F-statistic: 31.6 on 5 and 69 DF, p-value: < 2.2e-16Comprendre son modèle avec des graphiques
La fonction predict sert à calculer la valeur prédite
par le modèle pour chaque observation. Par exemple, ceci permet de
mettre en relation les valeurs observées et les valeurs prédites. La
fonction predict dans ce cas calcule la valeur prédite par
le modèle pour chaque observation dans la base de donnée.
plot(d$y, predict(m))
Soumettre de nouvelles valeurs à predict
Plus souvent, on veut illustrer l’effet de nos variables et pour
cela, il faut soumettre des valeurs à la fonction predict.
Cette fonction est très utile lorsque l’on veut un maximum de contrôle
sur les valeurs à prédire.
x <- seq(min(d$x1), max(d$x1), length.out = 50)
nd <- data.frame(x1 = x, x2 = mean(d$x2), x3 = mean(d$x3), fac = "B") # on fixe les autres variables à leur valeur moyenne
p <- predict(m, newdata = nd)
plot(d$x1, d$y) # observations
lines(x, p) # valeurs prédites
Package visreg
Le package visreg
permet d’illustrer facilement les prédictions d’un modèle (ou les effets
marginaux des différentes variables explicatives). Il utilise en
arrière-plan la fonction predict.
library(visreg)
par(mfrow = c(2, 2))
visreg(m)
Package ggeffects
Le package ggeffects
permet également d’illustrer les prédictions d’un modèle, mais ce
package se base sur l’utilisation du package ggplot2 pour construire les
graphiques. Le principe demeure le même: en arrière-plan la fonction
predict est utilisée. On peut également combiner les
différents graphiques générés en utilisant le package patchwork et
sa fonction wrap_plots.
library(ggeffects)
library(patchwork)
g <- ggpredict(m)
p <- plot(g, add.data = TRUE)
wrap_plots(p)
Interactions
Dans vos mots, qu’est-ce qu’une interaction?
Un exemple
# modèle modifié avec le facteur
m <- lm(y ~ x1 + x2 + x3 + x1 * fac, data = d)
summary(m)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x1 * fac, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.0406 -1.9485 0.1316 2.1335 14.2374
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 98.50371 2.15978 45.608 < 2e-16 ***
## x1 0.14204 0.03151 4.508 2.7e-05 ***
## x2 1.15443 0.18388 6.278 2.9e-08 ***
## x3 -2.44076 1.86546 -1.308 0.19521
## facB 0.66407 2.63788 0.252 0.80201
## facC -1.01455 3.01397 -0.337 0.73746
## x1:facB 0.01671 0.04332 0.386 0.70090
## x1:facC 0.13612 0.04871 2.794 0.00678 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.692 on 67 degrees of freedom
## Multiple R-squared: 0.7329, Adjusted R-squared: 0.705
## F-statistic: 26.26 on 7 and 67 DF, p-value: < 2.2e-16Illustration graphique
visreg(m, "x1", by = "fac", overlay = TRUE)
On a une interaction significative, mais est-elle importante d’un point de vue biologique? Peut-être pas… Toujours distinguer entre la signification statistique et la signification biologique.
Une valeur de p significative c’est bien beau, mais ce qui importe vraiment c’est la taille de l’effet.
Interactions avec ggeffects
g <- ggpredict(m, terms = c("x1", "fac"), message = FALSE)
plot(g, add.data = TRUE, jitter = FALSE)
Interprétion
Beaucoup plus facile avec un graphique qu’en regardant une table de coefficients.
En présence d’une interaction (significative), les effets simples sont plus difficilement interprétables et ne peuvent être interprétés sans faire référence à l’interaction (à moins d’ajustements particuliers voir Schielzeth 2010).
Si pour deux variables impliquées dans une interaction les effets simples ne sont pas significatifs, mais que leur interaction l’est, il faut considérer que ces deux variables ont un effet même si les coefficients associés aux effets simples ne sont pas significatifs.
Types d’erreurs
- Erreur de type I: rejeter l’hypothèse nulle, alors qu’elle est vraie (faux positif).
- Erreur de type II: ne pas rejeter l’hypothèse nulle, alors qu’elle est fausse (faux négatif).
Simulations
Simulons un modèle dans lequel les variables n’ont aucun effet.
Que se passera-t-il si on répète plusieurs fois ce scénario?
À quoi ressembleront nos valeurs de p?
Créons d’abord une fonction sim_model pour générer un
modèle linéaire fictif.
sim_model <- function(b0, b1, b2, b3) {
n <- 200
x1 <- runif(n, 0, 100)
x2 <- runif(n, 0, 10)
x3 <- runif(n, 0, 1)
mu <- b0 + b1 * x1 + b2 * x2 + b3 * x3
y <- rnorm(n, mu, 5)
d <- data.frame(y, x1, x2, x3)
m <- lm(y ~ x1 + x2 + x3, data = d)
m
}Simulons
Ensuite, simulons nsims modèles fictifs et emmagasinons
ces modèles dans une liste nommée list_models.
nsims <- 500
list_models <- lapply(1:nsims, function(i) {
sim_model(b0 = 100, b1 = 0, b2 = 0, b3 = 0)
})Si on illustre les valeurs de p obtenues pour chacun des coefficients associés aux variables \(x_{n}\) à l’aide d’un histogramme, à quoi ressemblera cet histogramme?
Valeurs de p obtenus
p <- lapply(list_models, function(i) {
summary(i)$coef[2:4, 4]
})
hist(unlist(p), breaks = seq(0, 1, by = 0.05), xlab = "Valeurs de p")
Conclusions
Dans une certaine proportion des cas ( ~ 5% ), on conclut à un effet, alors qu’il n’y en a pas pas! C’est le fameux seuil \(\alpha\) de 0.05.
Sous l’hypothèse nulle, les valeurs de p sont distribuées uniformément entre 0 et 1.
Modèles linéaires standards
\[y \sim \mathcal{N}(\mu, \sigma^2)\]
\[\mu=\beta_0 + \beta_1 X_1 + \beta_2 X_2 + ... + \beta_k X_k\]
Formulation matricielle
\[y \sim \mathcal{N}(\beta\mathbf{X}, \mathbf{I}\sigma^2)\]
\(\beta\) est un vecteur de coefficients.
\(\mathbf{X}\) est une matrice avec les valeurs des variables explicatives.
\(\mathbf{I}\) est une matrice identité qui a pour effet d’assigner une variance égale \(\sigma^2\) à toutes les observations.
Modèles linéaires généraux
\[y \sim \mathcal{N}(\beta\mathbf{X}, \mathbf{\Sigma})\]
\(\mathbf{\Sigma}\) est une matrice de variance-covariance qui permet beaucoup plus de flexibilité en terme de variance et de dépendance entre les observations. En général, des paramètres seront associés à la matrice \(\mathbf{\Sigma}\) pour permettre différentes situations. Dans le cas d’un modèle linéaire standard, on a \(\mathbf{\Sigma} = \mathbf{I}\sigma^2\).
Ex.: variance différente pour les niveaux d’un facteur, variance augmentant avec certaines valeurs de \(x\), autocorrélation temporelle ou spatiale, etc.
Ces différentes situations peuvent être obtenues entre autres avec
les packages nlme
(fonctions gls ou lme) ou glmmTMB.
Modèles linéaires généralisés (GLM)
\[y \sim \mathcal{ExpFam}(\theta, ...)\]
\[\mathcal{g}(\mu)=\beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k\]
Structure d’un GLM
\[y \sim \mathcal{ExpFam}(\theta, ...)\]
\[ \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k\]
\[ g(E(y)) = g(\mu) = \eta \]
Composante systématique
\[ \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k\]
Le prédicteur linéaire \(\eta\) permet de relier les variables explicatives à la variable réponse.
Composante aléatoire
\[ y \sim \mathcal{ExpFam}(\theta) \]
Les observations \(y\) suivent une distribution particulière appartenant à la famille des distributions exponentielles avec des paramètres \(\theta\). La distributions normale, de Poisson, Binomiale et Gamma font partie de la famille des distributions exponentielles.
Fonction de lien
\[ g(E(y)) = g(\mu) = \eta \]
La fonction de lien (link function) est une transformation qui permet de relier la réponse attendue \(E(y)\) (i.e. la moyenne) au prédicteur linéaire.
C’est ce qui permet de relier les variables explicatives à la variable réponse. Pour chaque distribution qui pourrait être utilisée dans le cadre d’un GLM, une fonction de lien dite canonique est utilisée par défaut, mais il existe plusieurs fonctions de liens possibles.
GLM Poisson
\[ y \sim \mathcal{Pois}(\lambda) \]
\[ log(\lambda) = \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k \]
On a donc des observations \(y\) qui suivent une distribution de Poisson avec un paramètre \(\lambda\). Une fonction de lien \(log\) permet de relier linéairement les variables explicatives à la valeur attendue \(E(y) = \lambda\) de la variable réponse.
Distribution de Poisson
Cette distribution est souvent utilisée pour décrire le nombre de fois qu’un événement survient dans un intervalle de temps ou d’espace.
\[f(k,\lambda) = Pr(y = k) = \frac{\lambda^{k}e^{-\lambda}}{k!}\]
Cette fonction permet de calculer la probabilité d’observer un événement \(k\) fois avec une distribution de Poisson ayant un paramètre \(\lambda\). La moyenne de cette distribution est égale à \(\lambda\).
Cette distribution a la propriété que \(\lambda = E(y) = Var(y)\). Autrement dit, la moyenne et la variance sont égales et correspondent au paramètre \(\lambda\) de la distribution.
D’où provient-elle?
Simulons une distribution de Poisson en distribuant aléatoirement 200 observations dans un carré 10 x 10. Ensuite, divisons ce carré en 100 cellules.
library(sf)
n <- 200
x <- runif(n, 0, 10)
y <- runif(n, 0, 10)
sfc <- st_sfc(st_polygon(list(rbind(c(0, 0), c(10, 0), c(10, 10), c(0, 0)))))
grid <- st_make_grid(sfc, cellsize = 1)
plot(grid, axes = TRUE)
points(x, y)
Nb d’observations par cellule
Ensuite, calculons le nombre d’observations par cellule et illustrons la distribution du nombre d’observations pour chaque cellule. En moyenne, on devrait retrouver 2 observations par cellule (200 obs. / 100 cell. = 2).
o <- st_intersects(grid, st_as_sf(data.frame(x, y), coords = c("x", "y")))
co <- sapply(o, length)
hist(co, breaks = seq(-0.5, max(co) + 0.5, by = 1), freq = FALSE, main = "", xlab = "Nb d'observations")
lines(0:max(co), dpois(0:max(co), lambda = 2), type = "b", cex = 2)
legend("topright", lty = c(NA, 1), legend = c("Observée", "Théorique"), bty = "n", cex = 1, pch = c(22, 21), pt.bg = c("gray", "white"))
Ce serait la même chose si on simulait des observations aléatoires dans le temps et si on comptait le nb d’observations par intervalles de temps.
Simulons un modèle
Simulons un modèle à partir d’une distribution de Poisson.
n <- 100
x <- runif(n, 0, 10)
lambda <- exp(-2 + 0.4 * x)
y <- rpois(n, lambda)
d <- data.frame(y, x)
m <- glm(y ~ x, family = "poisson", data = d)
summary(m)##
## Call:
## glm(formula = y ~ x, family = "poisson", data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.46287 0.30996 -7.946 1.93e-15 ***
## x 0.45294 0.03784 11.971 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 289.290 on 99 degrees of freedom
## Residual deviance: 85.894 on 98 degrees of freedom
## AIC: 245.83
##
## Number of Fisher Scoring iterations: 5Coefficients
Les coefficients estimés sont très similaires à ce qui a été utilisé pour générer le modèle (-2 et 0.4).
coef(m)## (Intercept) x
## -2.4628709 0.4529379Prédictions
Illustrons l’effet de la variable explicative.
visreg(m)
Sous quelle échelle?
visreg(m, scale = "response")
Échelle du prédicteur linéaire ou de la réponse
par(mfrow = c(1, 2))
visreg(m)
visreg(m, scale = "response")
Remarquez que dans le graphique de gauche, les valeurs en y sont négatives ou positives. En fait, la fonction log permet de transformer les comptes allant de 0 à l’\(+\infty\) en valeurs allant de \(-\infty\) à \(+\infty\). Sous l’échelle log, les valeurs entre 0 et 1 sont négatives et les valeurs supérieures à 1 sont positives.
\[ log(1) = 0 \] \[ log(0.1) = -2.302586 \] \[ log(5) = 1.609438 \]
ggeffects
g <- ggpredict(m, terms = c("x [n=50]"))
plot(g, add.data = TRUE, jitter = FALSE)
predict
v <- seq(min(x), max(x), length.out = 50)
plot(y ~ x, data = d)
p <- predict(m, data.frame(x = v))
lines(v, p)
Argument type
Par défaut, l’argument type = "link" et il faut
spécifier type = "response" pour obtenir les prédictions
sous l’échelle de la réponse.
v <- seq(min(x), max(x), length.out = 50)
plot(y ~ x, data = d)
p <- predict(m, data.frame(x = v), type = "response")
lines(v, p)
Suppositions de base
On pourrait être tenté d’utiliser la fonction plot pour
vérifier les conditions d’applications.
par(mfrow = c(2, 2), mar = c(4, 4, 2, 2))
plot(m)
Quelles sont les suppositions de bases?
\[ y \sim \mathcal{Pois}(\lambda) \]
\[ log(\lambda) = \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k \]
Certaines suppositions sont les mêmes que pour un modèle linéaire standard avec une distribution, mais d’autres diffèrent.
Vérification avec DHARMa
Une façon de vérifier les supposition de bases des modèles est d’utliser le modèle pour simuler des observations et ensuite comparer ces observations simulées aux observations réelles. Si les observations réelles ne ressemblent pas aux observations simulées à partir du modèle, c’est que le modèle ne représente pas bien ce qui se passe avec les données réelles et donc que qqch devrait être modifié dans le modèle.
Le package DHARMa se base sur ces simulations pour comparer les observations attendues par le modèle aux observations réelles. Plus spécifiquement, la méthode utilise les résidus quantiles pour comparer chaque observation aux observations attendues pour cette observations. Les résidus quantiles sont rapportés entre 0 et 1 et si le modèle est approprié pour les données, on devrait s’attendre à ce que les résidus soient distributés uniformément entre 0 et 1 pour toutes les observations.
library(DHARMa)
s <- simulateResiduals(m)
plot(s)
Si le modèle est appropré pour les données, on devrait voir ici une distribution uniforme des points et il ne devrait pas y avoir de patrons particuliers dans la distribution des points dans le graphique de droite. En particulier, les lignes de prédictions illustrées devraient être alignées sur les trois lignes représentant les quantiles 25, 50 et 75%. Dans le graphique quantile-quantile de gauche, les points devraient être alignés sur la droite.
Autres vérifications
Les problèmes de sous- ou sur- dispersion et de surabondance de 0 (zéro-inflation) sont très fréquents avec les GLM. Ces deux fonctions de DHARMa permettent d’évaluer si les observations réelles diffèrent des observations simulées en terme de présence de 0 ou de dispersion. En rouge, ce sont les observations réelles et en gris les observations simulées. Si les observations réelles sont bien en dehors des observations simulées, c’est qu’il y a un problème avec le modèle et il faut un ajustment.
par(mfrow = c(1, 2))
testDispersion(s)
testZeroInflation(s)
GLM binomial
\[ y \sim \mathcal{Binom}(n,p) \]
\[ logit(p) = log\left(\frac{p}{1-p}\right) = \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k \]
Régression logistique
\[ y \sim \mathcal{Binom}(1,p) = \mathcal{Bernoulli}(p)\]
La distribution de Bernoulli est un cas particulier de la distribution binomiale avec un tirage (e.g. un individu est vivant ou mort, un individu est mâle ou femelle, etc.)
Distribution binomiale
\[f(k,n,p) = Pr(y = k) = \binom{n}{k}p^k(1-p)^{n-k}\]
où
\[\binom{n}{k}=\frac{n!}{k!(n-k)!}\]
Ce qui représente le nombre de façons de prendre \(k\) éléments parmi \(n\) éléments.
Cette fonction permet de calculer la probabilité d’observer \(k\) succès si \(n\) tirages sont faits et que la probabilité de succès est de \(p\)
Cette distribution a la propriété que \(E(y) = np\) et \(Var(y) = np(1-p)\) , autrement dit, la moyenne et la variance sont entièrement déterminées par la probabilité de succès et par le nombre de tirages. On a donc une autre distribution assez restrictive.
Exemple
Supposons qu’on tire pile (0) ou face (1) 10 fois et qu’on calcule le nombre de face et qu’on répète cette exercice 500 fois. Quelle sera la distribution du nombre de faces obtenus?
n <- 500
co <- sapply(1:n, function(i) {
s <- sample(0:1, 10, prob = c(0.5, 0.5), replace = TRUE)
sum(s) # on additione les faces (1)
})
hist(co, breaks = seq(-0.5, max(co) + 0.5, by = 1), freq = FALSE, main = "", xlab = "Nb de faces")
lines(0:max(co), dbinom(0:max(co), prob = 0.5, size = 10), type = "b", cex = 2)
legend("topright", lty = c(NA, 1), legend = c("Observée", "Théorique"), bty = "n", cex = 1, pch = c(22, 21), pt.bg = c("gray", "white"))
Simulons un modèle
logit <- function(p) {
log(p/(1 - p))
}
inv.logit <- function(x) {
exp(x)/(1 + exp(x))
}
n <- 200
x <- runif(n, 0, 10)
z <- -2 + 0.5 * x
y <- rbinom(n, size = 1, prob = inv.logit(z))
d <- data.frame(y, x)
m <- glm(y ~ x, family = "binomial", data = d)
summary(m)##
## Call:
## glm(formula = y ~ x, family = "binomial", data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.68335 0.43063 -6.231 4.63e-10 ***
## x 0.63864 0.08853 7.214 5.45e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 272.74 on 199 degrees of freedom
## Residual deviance: 185.42 on 198 degrees of freedom
## AIC: 189.42
##
## Number of Fisher Scoring iterations: 5Visualisation
Avec visreg
par(mfrow = c(1, 2))
visreg(m)
visreg(m, scale = "response")
Remarquez que dans le graphique de gauche, les valeurs en y sont négatives ou positives. En fait, la fonction logit permet de transformer les probabilités allant de 0 à 1 en valeurs allant de \(-\infty\) à \(+\infty\). Sous l’échelle logit, une valeur de 0 corrrespond à une probabilité de 0.5.
\[ log\left(\frac{p}{1-p}\right) = log\left(\frac{0.5}{1-0.5}\right) = log(1) = 0 \]
Visualisation
Avec ggeffects
g <- ggpredict(m, terms = c("x [n=50]"))
plot(g, add.data = TRUE)
Suppositions de base
par(mfrow = c(2, 2))
s <- simulateResiduals(m)
testUniformity(s)
testQuantiles(s)
testDispersion(s)
testZeroInflation(s)
Problèmes
Les problèmes les plus fréquents sont la surdispersion et la surabondance de zéros (zero-inflation)
Surdispersion
C’est lorsque les données sont plus variables que ce qui est attendu par le modèle. Les distribution de Poisson et binomiale sont très restrictives par rapport à leur variance et les problèmes de surdispersion sont très fréquents.
En général, une surdispersion ignorée fera en sorte que l’incertitude sera sous-estimée (trop de confiance dans nos résultats).
Plusieurs causes possibles (variables explicatives manquantes, phénomènes d’aggrégation, hétérogénéité, etc.)
Solutions:
Poisson -> Négative binomiale (glm.nb,
glmmTMB, brms, etc.)
Binomiale -> Béta binomiale (glmmTMB,
brms, etc.)
Obervation-Level Random Effects (OLRE) (voir Harrison 2015 et Harrison 2014).
Exemple
Voici un exemple de GLM de Poisson avec de la surdispersion. Pour créer cette surdispersion, j’utilise une distribution négative binomiale pour générer les observations au lieu d’une distribution de Poisson.
n <- 100
x <- runif(n, 0, 10)
lambda <- exp(-2 + 0.4 * x)
y <- rnbinom(n, mu = lambda, size = 0.5) #
d <- data.frame(y, x)
m <- glm(y ~ x, family = "poisson", data = d)
g <- ggpredict(m, terms = "x")
plot(g, add.data = TRUE)
Vérification avec DHARMa
Les résidus ne sont clairement pas distributés uniformément sur le second graphique et la surdispersion est évidente sur le 3e graphique.
s <- simulateResiduals(m)
par(mfrow = c(1, 4))
testUniformity(s)
testQuantiles(s)
testDispersion(s)
testZeroInflation(s)
Correction avec glmmTMB
library(glmmTMB)
m <- glmmTMB(y ~ x, family = "nbinom1", data = d)
s <- simulateResiduals(m)
par(mfrow = c(1, 4))
testUniformity(s)
testQuantiles(s)
testDispersion(s)
testZeroInflation(s)
Sousdisperion
C’est lorsque les données sont moins variables que ce qui est attendu par le modèle. Ceci est beaucoup plus rare et probablement moins dommageable (e.g. taille de couvée/portée chez les oiseaux/mammifères voir Brooks et al. 2019).
En général, une sous-dispersion ignorée fera en sorte que l’incertitude sera surestimée (trop conservateurs dans nos résultats) ce qui est un moindre mal.
Solutions:
Poisson -> glmmTMB, family = "compois"
Zéro-inflation
Surabondance de 0 dans les observations par rapport à ce qui est attendu par le modèle. Il est possible aussi d’avoir une sous-abondance de 0 et/ou encore une surabondance de 1, mais ce dernier cas est un peu plus rare.
Solutions:
modèle zero-inflated
(glmmTMB,pscl,brms, etc.)
À noter qu’une surdipersion peut être due à une surabondance de 0 et vice-versa. Ce sont deux problèmes qui peuvent interagir.
Extensions
Gamma \(]0,\infty\)
Variable réponse continue, strictement positive et avec une asymétrie à droite.
Souvent une alternative à une transformation \(log\).
Ex.: concentration d’un composé, biomasse, etc.
Tweedie \([0,\infty\)
Variable réponse continue positive et incluant des 0.
Ex.: biomasse capturée, précipitations, etc.
Voir Lecomte et al. 2013
Beta \(]0,1[\) ou \([0,1]\)
Variable continue entre 0 et 1 exclusivement.
Pour proportion ou % ne provenant pas d’un ratio entre entiers.
Si 0 et 1 possibles, certaines transformations peuvent être utilisées.
Ex.: proportion d’une surface affectée, proportion de temps passé en alimentation, etc.
Voir Douma & Weedon 2019
2. Modèles mixtes/hiérarchiques
\[ y = \alpha + \beta x +
\epsilon \]
ou
\[ \mu = \alpha + \beta x
\]
\[ y \sim \text{Normale}(\mu,\sigma^{2}) \]
Modèle linéaire hiérarchique simple
\[ y_{ij} = \alpha + \beta x_{ij}
+ \tau_{i} + \epsilon_{ij}\]
\[ \tau_{i} \sim \text{Normale}(0,\,\sigma_{\tau}^{2}) \]
\[ \epsilon_{ij} \sim \text{Normale}(0,\sigma_{\epsilon}^{2})\ \]
où il y a \(i\) groupes et \(j\) observations par groupe
Autre formulation
\[ \mu_{ij} = \alpha_{i} + \beta
x_{ij}\]
\[ \alpha_{i} \sim \text{Normale}(\alpha,\,\sigma_{\alpha}^{2})\ \]
\[ y_{ij} \sim \text{Normale}(\mu_{ij},\,\sigma_{\epsilon}^{2})\ \]Un modèle est dit hiérarchique lorsque certains de ses paramètres sont définis par une distribution de probabilité.
Simulons un exemple
set.seed(12345)
ng <- 20
nobs <- 10
n <- ng * nobs
id <- gl(ng, nobs)
# valeurs fictives des variables explicatives selon une distribution uniforme
x <- runif(n, 0, 100)
# valeurs des paramètres
a <- 100
b <- 1
# on génère l'effet aléatoire pour chaque groupe
rea <- rep(rnorm(ng, 0, 20), each = nobs)
# prédicteur linéaire
mu <- (a + rea) + b * x
# observations à partir d'une distribution normale avec un écart-type de 5
y <- rnorm(n, mu, 5)
# données
d <- data.frame(y, x)Illustrons le modèle
library(lme4)
library(ggeffects)
m <- lmer(y ~ x + (1 | id), data = d)
g <- ggpredict(m, terms = c("x", "id [sample=20]"), type = "random")
plot(g, add.data = TRUE, show_ci = FALSE, colors = rep("black", 20), show_legend = FALSE)
Illustrons les pentes par individus
g <- ggpredict(m, terms = c("x", "id [sample=20]"), type = "random")
plot(g, add.data = TRUE, show_ci = FALSE, colors = sample(colors(), 20))
Retrouvons les paramètres
La variance de l’effet aléatoire était de 20² et la variance résiduelle était de 5².
summary(m)## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x + (1 | id)
## Data: d
##
## REML criterion at convergence: 1314.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.50905 -0.67270 -0.04039 0.61012 2.58655
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 386.9 19.67
## Residual 25.3 5.03
## Number of obs: 200, groups: id, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 101.35098 4.46512 22.70
## x 0.99676 0.01262 78.97
##
## Correlation of Fixed Effects:
## (Intr)
## x -0.153Récupérons ces effets
fixef(m)## (Intercept) x
## 101.3509752 0.9967622ranef(m)## $id
## (Intercept)
## 1 2.8823702
## 2 -22.7874952
## 3 5.4199340
## 4 -28.1373544
## 5 1.0787346
## 6 -10.6043879
## 7 -6.4207965
## 8 31.3275567
## 9 -10.5628148
## 10 4.4707573
## 11 -24.5843971
## 12 -26.8974001
## 13 24.9063468
## 14 23.4502635
## 15 -2.0989706
## 16 -12.2974161
## 17 -1.4157682
## 18 8.9803887
## 19 44.1496752
## 20 -0.8592262
##
## with conditional variances for "id"apply(ranef(m)$id, 2, sd)## (Intercept)
## 19.60505Coefficients
coef(m)## $id
## (Intercept) x
## 1 104.23335 0.9967622
## 2 78.56348 0.9967622
## 3 106.77091 0.9967622
## 4 73.21362 0.9967622
## 5 102.42971 0.9967622
## 6 90.74659 0.9967622
## 7 94.93018 0.9967622
## 8 132.67853 0.9967622
## 9 90.78816 0.9967622
## 10 105.82173 0.9967622
## 11 76.76658 0.9967622
## 12 74.45358 0.9967622
## 13 126.25732 0.9967622
## 14 124.80124 0.9967622
## 15 99.25200 0.9967622
## 16 89.05356 0.9967622
## 17 99.93521 0.9967622
## 18 110.33136 0.9967622
## 19 145.50065 0.9967622
## 20 100.49175 0.9967622
##
## attr(,"class")
## [1] "coef.mer"Complexifions le modèle
b## [1] 1reb <- rep(rnorm(ng, 0, 0.5), each = nobs)
mu <- (a + rea) + (b + reb) * x
y <- rnorm(n, mu, 5)
d <- data.frame(y, x)Illustrons les pentes
m <- lmer(y ~ x + (x | id), data = d)
g <- ggpredict(m, terms = c("x", "id [sample=20]"), type = "random")
plot(g, add.data = TRUE, show_ci = FALSE, colors = rep("black", 20), show_legend = FALSE)
Retrouvons les paramètres
La variance de l’effet aléatoire sur la pente était de (0.5)².
summary(m)## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x + (x | id)
## Data: d
##
## REML criterion at convergence: 1370.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.19310 -0.58325 0.02874 0.58522 2.42595
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 388.2872 19.7050
## x 0.2105 0.4588 -0.27
## Residual 22.0018 4.6906
## Number of obs: 200, groups: id, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 100.4564 4.4699 22.474
## x 0.8959 0.1033 8.672
##
## Correlation of Fixed Effects:
## (Intr)
## x -0.283
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00784532 (tol = 0.002, component 1)Récupérons ces effets
fixef(m)## (Intercept) x
## 100.4563555 0.8958772ranef(m)## $id
## (Intercept) x
## 1 -4.259743 -0.11207853
## 2 -24.555894 0.09388618
## 3 9.760813 0.08930651
## 4 -23.972706 -0.43249104
## 5 6.379230 0.06853379
## 6 -10.672292 0.92265805
## 7 -4.656865 -0.14388477
## 8 33.097624 -0.22946144
## 9 -7.593899 0.23044679
## 10 -2.310431 0.38486838
## 11 -32.815439 -0.47826203
## 12 -26.583571 0.53612804
## 13 22.768061 0.34741573
## 14 22.280452 -0.58389607
## 15 1.432798 0.75717249
## 16 -8.201353 -0.31288440
## 17 1.063482 0.10990165
## 18 10.101182 0.03586551
## 19 39.745709 -0.88012487
## 20 -1.007156 -0.40309995
##
## with conditional variances for "id"apply(ranef(m)$id, 2, sd)## (Intercept) x
## 19.4359900 0.4557504coef(m)## $id
## (Intercept) x
## 1 96.19661 0.78379868
## 2 75.90046 0.98976339
## 3 110.21717 0.98518372
## 4 76.48365 0.46338617
## 5 106.83559 0.96441100
## 6 89.78406 1.81853525
## 7 95.79949 0.75199243
## 8 133.55398 0.66641577
## 9 92.86246 1.12632400
## 10 98.14592 1.28074559
## 11 67.64092 0.41761518
## 12 73.87278 1.43200525
## 13 123.22442 1.24329294
## 14 122.73681 0.31198114
## 15 101.88915 1.65304970
## 16 92.25500 0.58299281
## 17 101.51984 1.00577886
## 18 110.55754 0.93174272
## 19 140.20206 0.01575234
## 20 99.44920 0.49277726
##
## attr(,"class")
## [1] "coef.mer"Formulations
m <- lmer(y ~ x + (1 | id), data = d)
m <- lmer(y ~ x + (x | id), data = d) # équivalent à (1 + x|id)
m <- lmer(y ~ x + (0 + x | id), data = d)
m <- lmer(y ~ x + (1 | id) + (0 + x | id), data = d)Corrélations intercepts/pentes
m <- lmer(y ~ x + (1 | id) + (0 + x | id), data = d)
Le dernier est équivalent à
m <- lmer(y ~ x + (x || id), data = d)
et signifie qu’on assume que l’intercept et la pente ne
sont pas corrélées ou qu’on ne veut pas estimer cette corrélation. En
pratique, il n’y a pas énormément de raison de faire cela et par défaut
le modèle va estimer cette corrélation.
Effets aléatoires croisés
Supposons que nous avons des individus observés plusieurs fois par année et ce sur plusieurs années.
m <- lmer(y ~ x + (1 | id) + (1 | year), data = d)Effets aléatoires nichés
Supposons maintenant que les individus sont annuels et qu’ils ne sont
observés que pour des années particulières.
m <- lmer(y ~ x + (1 | year/id), data = d)
Équivalent à:
m <- lmer(y ~ x + (1 | year) + (1 | year:id), data = d)Quand utiliser les effets aléatoires?
Mesures répétées sur certaines unités ou individus
Les observations sont regroupées sous des structures ou ne sont pas indépendantes
Il y a de la réplication à plusieurs niveaux de la hiérarchie
On veut estimer la variance ou généraliser à une population d’unités d’échantillonnage
Effet spatial, temporel, phylogénétique
Autocorrélation spatiale, temporelle, phylogénétique, etc…
Régularisation (shrinkage)
Complete pooling
No pooling
Partial pooling
Un exemple
Effet du n
3. Vraisemblance et statistiques bayésiennes
Rappel des probabilités
Quelle est la probabilité d’obtenir 2 fois pile en tirant une pièce de monnaie 2 fois?
Pile ou face
On doit multiplier la probabilité de chaque événement pour obtenir la probabilité combinée de ces deux événements. Si \(P(pile) = 0.5\) (i.e. la probabilité d’obtenir pile), on a:
\[P(pile,pile) = P(pile) \times
P(pile) = 0.5 \times 0.5 = 0.25\]
Avec R
On peut calculer ceci avec R et la fonction dbinom.
dbinom(0, size = 1, p = 0.5) * dbinom(0, size = 1, p = 0.5)## [1] 0.25
On peut également appliquer la fonction dbinom à un
vecteur de valeur. Ceci nous donne un vecteur de (densité de)
probabilité et on peut multiplier ces valeurs entre elles avec la
fonction prod.
prod(dbinom(c(0, 0), size = 1, p = 0.5))## [1] 0.25Vraisemblance
La vraisemblance (likelihood) représente la plausibilité de
valeurs de paramètres en fonction de données observées. L’utilisation de
la vraisemblance implique l’utilisation de distribution de probabilités
pour représenter le lien entre les observations et les modèles.
Pour quantifier la vraisemblance, on peut utiliser une fonction de
vraisemblance (likelihood function).
\[L(\theta|data) = P(data|\theta)\]
\(data\): données
\(\theta\): paramètres du modèle
\(L(\theta|data)\): fonction de
vraisemblance
\(P(data|\theta)\):
vraisemblance (likelihood)
Ici, \(L(\theta|data)\) est la
fonction de vraisemblance et l’équation précédente nous dit que la
vraisemblance des paramètres conditonnelle aux paramètres \(\theta\) est égale à la (densité de)
probabilité des données conditionnelle aux paramètres \(\theta\).
Maximum de vraisemblance
L’idée derrière le maximum de vraisemblance (maximum likelihood) est de déterminer la valeur des paramètres qui rend l’observation des données la plus probable.
Une grande partie des paramètres que nous estimons sont
estimés par la méthode du maximum de vraisemblance.
Certaines
méthodes plus classiques comme la minimisation de la somme des carrés
dans la régression standard sont équivalentes à la maximisation de la
vraisemblance. La sélection de modèles par AIC repose sur l’utilisation
du maximum de vraisemblance.
Exemple
Simulons n = 20 présences/absences avec une probabilité
de p = 0.5. Ensuite, calculons la vraisemblance d’observer
les données générées en fonction de plusieurs valeurs du paramètre p,
qui ici est une probabilité.
set.seed(123)
p <- 0.5
y <- rbinom(20, 1, p) # tirage aléatoire de 20 valeurs de 0/1 avec une probabilité de 0.5
probs <- seq(0, 1, by = 0.01) # valeurs du paramètre qui seront évaluées
l <- sapply(probs, function(i) {
prod(dbinom(y, size = 1, prob = i)) # likelihood
})
Pour calculer la vraisemblance, on utilise la fonction dbinom
qui permet de calculer pour chaque valeur observée une (densité de)
probabilité. En calculant le produit prod de toutes ces
valeurs, on obtient la valeur de la fonction de vraisemblance pour
l’observation de l’ensemble de ces données pour une valeur donnée du
paramètre p. On calcule ensuite ces valeurs de
vraisemblance pour un ensemble de valeur de p.
Profil de vraisemblance
Par la suite, on peut illustrer la vraisemblance d’obtenir les
valeurs observées pour les différentes valeurs de p. Ceci
permet d’illustrer le profile de vraisemblance (likelihood
profile).
plot(probs, l, type = "l", ylab = "Likelihood", xlab = "Valeur du paramètre p")
abline(v = sum(y)/length(y), lty = 3)
abline(v = 0.5)
legend("topleft", lty = c(1, 3), legend = c("Probabilité réelle", "Proportion de 1 dans l'échantillon"), bty = "n")
probs[which.max(l)]## [1] 0.55
Le maximum de cette courbe est le maximum de vraisemblance
(maximum likelihood), i.e. la valeur de p pour
laquelle la probabilité d’observer les données est la plus grande. Dans
ce cas précis, cette valeur est la proportion de 1 dans les données. La
plupart des paramètres que nous estimons sont obtenus avec cette
méthode.
Avec un glm
Voici ce qu’on obtient avec un glm si on tente d’estimer cette probabilité. En principe, l’intercept représente cette probabilité, mais elle est sur l’échelle logit \(log(p/(1-p))\). Il faut donc retransformer le tout sous l’échelle d’une probabilité avec une fonction logit inverse.
m <- glm(y ~ 1, family = binomial)
summary(m)##
## Call:
## glm(formula = y ~ 1, family = binomial)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2007 0.4495 0.446 0.655
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27.526 on 19 degrees of freedom
## Residual deviance: 27.526 on 19 degrees of freedom
## AIC: 29.526
##
## Number of Fisher Scoring iterations: 3exp(coef(m))/(1 + exp(coef(m)))## (Intercept)
## 0.55Intervalle de confiance
On peut également estimer un intervalle de confiance autour de cette
probabilité avec la fonction confint. Différentes méthodes
existent pour obtenir un intervalle de confiance à partir d’un profil de
vraisemblance, mais nous n’arborderons pas cela ici (i.e.
profile likelihood confidence intervals).
ci <- confint(m) # calcul de l'intervalle de confiance sur les coefficients## Waiting for profiling to be done...ci <- c(coef(m), ci) # intercept est l'estimation de la probabilité p
c(exp(ci)/(1 + exp(ci))) # intervalle de confiance## (Intercept) 2.5 % 97.5 %
## 0.5500000 0.3359294 0.7519795Log de la vraisemblance
Pour des raisons pratiques, on va souvent maximiser le log
de la vraisemblance (log-likelihood) pour éviter des problèmes
numériques associés aux très petites valeurs. On parle donc souvent de
log-likelihood (LL, logLik, etc.).
Par
exemple, on peut comparer la vraisemblance calculée par le glm à celle
que nous avons calculée “à la mitaine” avec la fonction
dbinom pour la valeur du paramètre qui maximise cette
vraisemblance (p = 0.55).
prod(dbinom(y, size = 1, prob = 0.55))## [1] 1.05415e-06exp(logLik(m))## 'log Lik.' 1.05415e-06 (df=1)Probabilités conditionnelles
Voyons d’abord ce qu’est une probabilité conditionnelle.
\[P(A|B)\]
Ici, on représente la
probabilité de \(A\) considérant que
\(B\) s’est produit. Par exemple,
quelle est la probabilité d’être infecté par la Covid-19 si on reçoit un
test positif.
\[P(\textrm{infecté}|\textrm{test
positif})\]
On pourrait aussi s’intéresser à la
probabilité d’avoir un test positif si on est effectivement infecté.
\[P(\textrm{test
positif}|\textrm{infecté})\]
Règle de Bayes
Il s’agit simplement d’une règle en théorie des probabilités qui
permet de relier des probabilités conditionnelles.
\[P(A|B) = \frac{P(B|A)P(A)}{P(B)}\]
\(P(A)\) veut dire la probabilité de
\(A\), indépendemment de \(B\).
\(P(A|B)\) veut dire la probabilité de \(A\) sachant \(B\). C’est une probabilité
conditionnelle.
Interprétation
Un exemple classique de cette règle est le suivant. Supposons que
\(P(A)\) désigne la probabilité d’être
infecté par la Covid-19 et que \(P(B)\)
désigne la probabilité de tester positif à un test de Covid-19. Ici, on
a aussi \(P(A')\) qui est la
probabilité de ne pas être infecté et \(P(B')\) qui est la probabilité d’un
test négatif.
Si on a:
\(P(A) =
1/500 = 0.002\): proportion de personnes infectées
\(P(A') = 1-1/500 = 0.998\): proportion de personnes qui ne sont pas infectées
\(P(B|A) = 0.98\): probabilité d’un test positif si une personne est infectée
\(P(B|A') = 0.05\): probabilité d’un test positif si une personne n’est pas infectée
\(P(B) = P(B|A)P(A) + P(B|A')P(A')\): il faut additionner les possibilités associées au cas où une personne est infectée et au cas où une personne n’est pas infectée.
Quelle est la probabilité d’être infecté si on reçoit un test positif
(en assumant qu’on se fait tester pour une raison autre que d’avoir des
symptômes) ?
\[P(A|B) = \frac{P(B|A)P(A)}{P(B|A)P(A)+P(B|A')P(A')} = \frac{0.98 \times 0.002}{(0.98 \times 0.002) + (0.05 \times 0.998)} \approx 0.0378\]
Dans ce cas fictif, la probabilité d’avoir réellement la
Covid-19 si le test est positif n’est que d’environ 3.78 %
Statistiques bayésiennes
Les statisiques bayésiennes se basent sur cette règle pour faire une
inteprétation probabiliste des paramètres à partir des
données.
\[P(\theta|data) =
\frac{P(data|\theta)P(\theta)}{P(data)}\]
\(data\): données
\(\theta\): paramètres du modèle
\(P(\theta|data)\): distribution
postérieure (posterior)
\(P(data|\theta)\): vraisemblance
(likelihood)
\(P(\theta)\): distribution a priori
(prior)
\(P(data)\):
distribution marginale des données (marginal likelihood)
Distribution postérieure (posterior)
\[P(\theta|data) =
\frac{P(data|\theta)P(\theta)}{P(data)}\]
La
distribution postérieure \(P(\theta|data)\) est ce que l’on cherche à
estimer. Il s’agit d’une distribution de probabilité qui décrit les
probabilités des différentes valeurs possibles du paramètre.
Contrairement aux analyses fréquentistes, on peut directement parler de
la probabilité qu’un paramètre soit entre telle et telle valeur.
Vraisemblance (likelihood)
\[P(\theta|data) =
\frac{P(data|\theta)P(\theta)}{P(data)}\]
Comme vu
précédemment, la vraisemblance \(P(data|\theta)\) représente la probabilité
des données conditionnelles à des valeurs de paramètres \(\theta\).
Distribution a priori (priors)
\[P(\theta|data) =
\frac{P(data|\theta)P(\theta)}{P(data)}\]
Les
distributions a priori \(P(\theta)\) représentent nos croyances ou
nos connaissances sur les valeurs des paramètres avant de faire les
analyses. On les représente à l’aide de distributions de probabilité.
Ces croyances représentent ce à quoi on s’attend en termes de valeurs en
se basant sur nos connaissances avant de faire notre expérience.
Distribution marginale des données
\[P(\theta|data) =
\frac{P(data|\theta)P(\theta)}{P(data)} =
\frac{P(data|\theta)P(\theta)}{\int
P(data|\theta)P(\theta)d\theta}\]
Cette quantité \(P(data)\) (marginal likelihood)
est beaucoup plus difficile à comprendre intuitivement. Elle représente
en quelque sorte la probabilité des données pour l’ensemble des valeurs
de paramètres possibles (d’où l’intégrale). C’est la distribution
marginale des données ou la vraisemblance marginale. Elle permet
également de faire en sorte que la distribution postérieure soit une
distribution de probabilité (i.e. aire sous la courbe = 1).
Il
s’agit d’une constante qui est généralement très difficile à calculer.
Des méthodes comme le MCMC permettent d’éviter de la calculer
directement dans l’obtention de la distribution postérieure.
En résumé
\[Posterior = \frac{Likelihood \times Prior}{Marginal\;Likelihood}\]
Puisque le dénominateur est une constante et que dans bien des cas, on peut éviter de le calculer, on peut plus simplement écrire:
\[Posterior \propto Likelihood \times Prior\]
En d’autres mots, la distribution postérieure est
proportionnelle à la vraisemblance et aux distributions a
priori
On pourrait aussi dire qu’avec une approche bayésienne, on
met à jour nos croyances ou nos connaissances à l’aide de données.
Avantages
Interprétation probabiliste plus facile (e.g. intervalle de confiance fréquentiste vs. bayésien)
Facilite l’estimation de modèles complexes (e.g. modèle hiérarchique)
Intégration plus facile de l’incertitude à toutes les étapes du processus
Facilite l’estimation de quantités dérivées (e.g. prédictions ou combinaisons de prédictions)
Le fait d’avoir à spécifier des priors nous force à comprendre ce que nos paramètres représentent
Façon de penser plus naturelle, philosophiquement plus satisfaisante?
Désavantages ?
Souvent plus complexe (mais pas toujours)
Modèles plus longs à estimer
Moins de solutions clés en main?
Nécessité de choisir et de justifier ses priors
Exemple
Simulons d’abord un modèle avec des paramètres connus.
set.seed(1234)
n <- 100
x <- runif(n, 0, 10)
y <- 10 + 2 * x + rnorm(n, 0, 10)
d <- data.frame(y, x)Quels priors utiliser?
Il est toujours préférable de visualiser les priors et de
réfléchir à ce qu’ils représentent.
par(mfrow = c(1, 3))
curve(dnorm(x, 0, 50), from = -200, to = 200, main = "Intercept normal(0,50)", ylab = "Density")
curve(dnorm(x, 0, 5), from = -20, to = 20, main = "Pente normal(0,5)", ylab = "Density")
curve(dgamma(x, 2, 0.1), from = 0, to = 100, main = "Écart-type gamma(2,0.1)", ylab = "Density")
Par exemple, pour la variance/écart-type du modèle, il faut
idéalement un prior avec des valeurs strictement positives, car
les valeurs négatives ne font pas de sens pour la variance. Une des
possibilités est d’utiliser la distribution gamma.
Modèle avec brms
La façon la plus “clé en main” de faire un modèle bayésien est
d’utiliser le package brms. Ce package
permet de spécifier des modèles d’une façon très similaire aux packages
fréquemment utilisés comme nlme, lme4, glmmTMB, etc. En arrière-plan,
brms utilise le programme STAN qui
est un environnement de programmation pour la modélisation bayésienne.
library(brms)
m <- brm(y ~ x,
data = d,
cores = 3,
chains = 3,
family = "gaussian",
prior = c(
set_prior("normal(0,5)", class = "b", coef="x"),
set_prior("normal(0,50)", class = "Intercept"),
set_prior("gamma(2,0.1)", class = "sigma")
)
)Visualisation
library(ggeffects)
g <- ggpredict(m, terms = c("x"))
plot(g, add.data = TRUE)
Distribution postérieure des prédictions
newdata <- data.frame(x = seq(min(d$x), max(d$x), length.out = 100))
pp <- posterior_epred(m, newdata, ndraws = 500)
plot(d$x, d$y, ylim = range(c(pp, d$y)), pch = 16, col = gray(0, 0.35))
lapply(1:nrow(pp), function(i){
lines(newdata$x, pp[i, ], col = adjustcolor("black", 0.05))
}) |> invisible()
Formulation du modèle
On a la composante aléatoire, la composante systématique et les priors.
\[y \sim
\mathcal{N}(\mu,\sigma^2)\] \[\mu =
\beta_{0}+\beta_{1}x\] \[\beta_{0}
\sim \mathcal{N}(0,50)\] \[\beta_{1}
\sim \mathcal{N}(0,5)\] \[\sigma \sim
\mathcal{Gamma}(2,0.1)\]
Sommaire
summary(m)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: y ~ x
## Data: d (Number of observations: 100)
## Draws: 3 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 3000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 10.97 1.85 7.25 14.52 1.00 2892 2159
## x 1.96 0.36 1.26 2.69 1.00 2997 2243
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 9.68 0.72 8.42 11.19 1.00 2611 1835
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Notez encore une fois ici qu’on peut faire une interprétation
probabiliste plus naturelle des intervalles de confiance. Contrairement
aux analyses fréquentistes, on peut dire que le probabilité que la
valeur d’un paramètre soit entre les bornes inférieures et supérieures
d’un intervalle de confiance à 95% est de 95%. En fait, avec une
approche bayésienne, on parle plutôt d’intervalle de crédibilité
(credible interval).
Visualisation des paramètres
plot(m)
Les distributions postérieures des paramètres sont illustrées à gauche. À droite, les chaînes de MCMC sont illustrées (nous y reviendrons).
MCMC
Le MCMC (Markov Chain Monte Carlo) est une classe
d’algorithmes permettant d’échantillonner à partir d’une distribution de
probabilité. Il s’agit d’une approche qui permet d’éviter des calculs
complexes, voire impossible, en utilisant une approche
d’échantillonnage.
Dans un contexte bayésien, le MCMC permet
d’échantillonner à partir d’une distribution postérieure. Autrement dit,
la distribution postérieure de chaque paramètre estimé dans un contexte
bayésien peut entre autre être obtenue par cette méthode
d’échantillonnage.
Sous certaines conditions, les valeurs
échantillonnées représenteront la distribution postérieure de chaque
paramètre.
Échantillonnage des chaînes
L’échantillonnage par MCMC se fait de façon séquentielle,
i.e. que les valeurs seront échantillonnées les unes à la suite
des autres. Dans bien des cas, les valeurs actuelles seront dépendentes
(corrélées) des valeurs précédentes. On peut parler ici de chaînes pour
désigner l’ensemble des valeurs échantillonnées de façon séquentielle.
Pour démarrer chaque chaîne, il faut d’abord sélectionner une
valeur de départ pour chaque paramètre estimé. Cette valeur de départ
peut potentiellement être très loin d’une valeur probable.
Éventuellement, les chaînes vont s’éloigner de ces valeurs de départ
pour converger vers les distributions postérieures. À ce moment, on dira
que l’échantillonnage se fait véritablement à partir des distributions
postérieures.
plot(m)
Convergence des chaînes
Avec les analyses bayésiennes où les paramètres sont estimés par
MCMC, il faut donc vérifier la convergence des chaînes pour s’assurer
que les valeurs échantillonnées représentent bien les distributions
postérieures.
On élimine souvent une certaine portion du début
des chaînes pour enlever les valeurs échantillonnées, alors que les
chaînes n’ont pas convergé. On parle souvent de burnin pour
désigner cette partie des chaînes.
Il est également très
fréquent d’utiliser plusieurs chaînes (3-4) pour un même paramètre ce
qui permet de s’assurer que les chaînes ont toutes convergé au même
endroit. En général, on choisira des valeurs de départ différentes pour
chacune des chaînes.
On peut vérifier la convergence des
chaînes avec entre autres le Rhat (potential scale
reduction statistic ou Gelman-Rubin statistic).
Intuitivement, pour un paramètre donné, cette valeur compare la variance
inter-chaîne à la variance intra-chaîne pour évaluer si les chaînes ont
convergé au même endroit. Si c’est le cas, le Rhat sera
très près de 1 (Rhat <= 1.01)
summary(m)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: y ~ x
## Data: d (Number of observations: 100)
## Draws: 3 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 3000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 10.97 1.85 7.25 14.52 1.00 2892 2159
## x 1.96 0.36 1.26 2.69 1.00 2997 2243
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 9.68 0.72 8.42 11.19 1.00 2611 1835
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Priors
L’utilisation de priors est souvent vue comme étant la
partie la plus controversée des statistiques bayésiennes. Qu’en
pensez-vous?
Idéalement, il faut toujours avoir une bonne idée
de l’influence des priors sur nos analyses. Pour cela il est
bien de comparer les distributions postérieures aux priors. On
peut également faire dans certains cas des analyses de sensibilité pour
quantifier cette influence.
Plus on a de données, moins les
priors auront d’influence sur l’estimation des paramètres.
Types de priors
Informatif: basé sur la théorie, reflète les connaissances.
Non-informatif: reflète l’absence de connaissances ou le désir de ne pas influencer les résultats.
Plat (flat): similaire à non-informatif.
Légèrement informatif (weaky informative): la théorie n’est pas assez développée pour utiliser des priors informatifs, mais assez pour savoir que certaines valeurs ne sont pas probables.
Par défaut: valeurs par défaut utilisées par les logiciels ou les packages.
Régularisant (regularizing): permet de circonscrire les valeurs obtenues et de pénaliser les valeurs élevées dans un contexte prédictif.
Etc.
Voir notamment Lemoine 2019 et Banner et al. 2020 pour une discussion sur l’utilisation des priors en écologie.
Quels priors utiliser?
Il faut aussi sélectionner une distribution approriée pour
représenter nos priors et les distributions utilisées doivent
être appropriées pour un paramètre donné. Le groupe qui travaille sur la
plateforme STAN a plusieurs
recommandations de priors à utiliser dépendemment du contexte
(priors pour des pentes, priors pour des variances,
priors pour des variances dans un contexte de modèles
hiérarchiques, etc.)
Voir Prior Choice Recommendations pour des recommandations.
Attention aux priors “non-informatifs”!
Les priors non-informatifs ne le sont pas toujours sous
l’échelle du paramètre. Par exemple, si vous tentez d’estimer une
probabilité et que vous utiliser une distribution normale comme
prior pour l’intercept du modèle, il se peut que ce
prior soit très informatif sous l’échelle du prédicteur
linéaire qui est relié à la probabilité par une fonction logit.
Supposons un glm binomial avec seulement un intercept.
\[ y \sim \mathcal{Binom}(1,p) =
\mathcal{Bernoulli}(p)\]
\[
logit(p) = log\left(\frac{p}{1-p}\right) = \beta_{0} \]
\[ \beta_{0} \sim \mathcal{N}(0,10) \]
Attention aux priors “non-informatifs”!
n <- 10000
x <- rnorm(n, 0, 10)
par(mfrow = c(1, 2))
hist(x, breaks = 50, main = "Valeurs a priori de l'intercept")
hist(exp(x)/(1 + exp(x)), breaks = 50, main = "Lorsque converties en probabilité")
Autrement dit, le prior est très informatif et pousse les
valeurs vers des probabilités extrêmes!
Étudier ses priors
La distribution prédictive préalable (Prior predictive distribution) permet de visualiser l’effet des priors sur l’échelle de la réponse.
m2 <- update(m, sample_prior = "only")newdata <- data.frame(x = seq(min(d$x), max(d$x), length.out = 100))
pp <- posterior_epred(m2, newdata, ndraws = 500)
plot(d$x, d$y, ylim = range(c(pp, d$y)), pch = 16, col = gray(0, 0.35))
lapply(1:nrow(pp), function(i){
lines(newdata$x, pp[i, ], col = adjustcolor("black", 0.05))
}) |> invisible()
Avec un prior aberrant
newdata <- data.frame(x = seq(min(d$x), max(d$x), length.out = 100))
m3 <- update(m,
sample_prior = "only",
prior = c(
set_prior("normal(50,20)", class = "b", coef="x"),
set_prior("normal(0,50)", class = "Intercept"),
set_prior("gamma(2,0.1)", class = "sigma")
),
)pp <- posterior_epred(m3, newdata, ndraws = 500)
plot(d$x, d$y, ylim = range(c(pp, d$y)), pch = 16, col = gray(0, 0.35))
lapply(1:nrow(pp), function(i){
lines(newdata$x, pp[i, ], col = adjustcolor("black", 0.05))
}) |> invisible()
Avec un modèle binomial
logit <- function(p) {
log(p / (1 - p))
}
inv.logit <- function(x) {
exp(x) / (1 + exp(x))
}
n <- 200
x <- runif(n, 0, 10)
z <- -2 + 0.5 * x
y <- rbinom(n, size = 1, prob = inv.logit(z))
d<-data.frame(y,x)
m <- brm(y ~ x,
data = d,
cores = 3,
chains = 3,
family = "bernoulli",
prior = c(
set_prior("normal(0,5)", class = "b", coef="x"),
set_prior("normal(0,5)", class = "Intercept")
)
)Distribution postérieure des prédictions
newdata <- data.frame(x = seq(min(d$x), max(d$x), length.out = 100))
pp <- posterior_epred(m, newdata, ndraws = 500)
plot(d$x, d$y, ylim = range(c(pp, d$y)), pch = 16, col = gray(0, 0.35))
lapply(1:nrow(pp), function(i){
lines(newdata$x, pp[i, ], col = adjustcolor("black", 0.05))
}) |> invisible()
Ce qu’indique notre prior
m2 <- update(m, sample_prior = "only")pp <- posterior_epred(m2, newdata, ndraws = 500)
plot(d$x, d$y, ylim = range(c(pp, d$y)), pch = 16, col = gray(0, 0.35))
lapply(1:nrow(pp), function(i){
lines(newdata$x, pp[i, ], col = adjustcolor("black", 0.05))
}) |> invisible()
Prédiction pour x = 2
pp <- posterior_epred(m2, data.frame(x = 2), ndraws = 500)
hist(pp, breaks = 50)
Un prior plus raisonnable?
m3 <- update(m,
sample_prior = "only",
prior = c(
set_prior("normal(0,1)", class = "b", coef="x"),
set_prior("normal(0,1)", class = "Intercept")
)
)pp <- posterior_epred(m3, newdata, ndraws = 500)
plot(d$x, d$y, ylim = range(c(pp, d$y)), pch = 16, col = gray(0, 0.35))
lapply(1:nrow(pp), function(i){
lines(newdata$x, pp[i, ], col = adjustcolor("black", 0.05))
}) |> invisible()
Recommandations sur les priors
Prior Choice Recommendations Recommandation par les développeurs de STAN
The use of Bayesian priors in Ecology: The good, the bad and the not great
Choosing priors in Bayesian ecological models by simulating from the prior predictive distribution
Outils
Il existe plusieurs outils pour faire des modèles bayésiens.
WinBUGS
Un des premiers outils, mais un peu lent
JAGS
Similaire à BUGS, mais un peu plus rapide
STAN Rapide et bien
documenté
Nimble Relativement
nouveau
MCMCglmm
Modèle animal?
brms
Facile d’utilisation et basé sur STAN
INLA Modèles
spatiaux